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This paper presents a simple and fast algorithm for computing the union, 
the intersection, the difference and the symmetrical difference of polygonal 
regions in the plane. The algorithm explicitly copes with degenerate cases and 
vertices of high degree so that the output of the algorithm satisfies its input 
restrictions. An expected running time of the algorithm is 
O(n log*n + k + z log n), where n (resp. z) is the total number of edges (resp. 
contours) in polygons describing input regions and k is the number of 
intersections between the edges. The presented algorithm outperforms most of 


known algorithms for polygon set operations and is relatively easy to 
implement. 


INTRODUCTION 


Boolean set operations on polygons are essential in various applications 
such as computer graphics, GIS, CAD/CAM, circuit design etc. The result of a 
set operation on even two simple convex polygons can possibly be a set of 
concave self-touching polygons with holes (see Fig. 1). 


Unfortunately, the traditional 
algorithms either do not cope with 
degenerate cases at all or break 
complex output polygons into simple 
ones, making the operations on the 
resulting polygon(s) impossible 
without additional preprocessing 
steps. Nevertheless, in various GIS 
and CAD applications there is a need 
for a closed set of polygon 
operations. This means that the output of an algorithm must satisfy its input 
restrictions. It is also desireable to allow holes and vertices of degree higher 
than two. Let us make a brief survey of the existing algorithms for polygon set 
operations. 


Figure 1 


Weiler and Atherton [5. , 12. , 13. ] proposed the following algorithm. 
First, all intersections of input polygon edges are determined. A complex set of 
rules is used for the different classes of intersections to rearrange the contours 
so that none of them intersects. Then each contour is traversed to collect the 
information about its owner(s). Finally, the nesting structure of the contours is 
determined and the resulting polygon subset is selected from this structure. The 
main disadvantages of this algorithm are: 


e the result of an operation must not necessarily consist of the 
minimum contours (see Table 4, A + B), 


e the problem of handling and describing self-touching polygons is not 
considered at all, though the algorithm can yield these polygons as its 
output, 


e the output of the algorithm can contain vertices of high degree, which 
are not allowed on its input, and it is extremely difficult to extend the 
algorithm to cope with such cases. 


Schutte [9. , 10. ] modified the Weiler algorithm by much more clear 
division into steps and the different algorithm for edge labeling. However, his 
algorithm has stronger input restrictions. 


The above algorithms give a non-closed set of polygon operations, 
because their output can contain self-touching contours (which are not allowed 
by both of them) and holes, which are not allowed by the Schutte algorithm. 


In this paper we describe a simple and fast algorithm which is free of the 
disadvantages mentioned above. The paper is organized as follows. In section 
2 the basic definitions and denotations are given, section 3 describes the 
algorithm itself, in sections 4 and 5 we give the analysis of the algorithm 
efficiency and discuss some implementation details. 


BASIC DEFINITIONS 


Let E(a, b) denote the edge starting at point a and ending at point b, 
which is not coincident with a. For b (resp. a) we will call the edge E previous 
(resp. next). A set of points x satisfying the conditions 0 < Z(bax) < ọ and 
2n-@ < Z(abx) < 27 is called the left neighborhood of edge E(a, b). Similarly, 
a set of points x satisfying the conditions 0< Z(abx)<@ and 
2n-@ < Z(bax) < 27 is called the right neighborhood of edge E(a, b) (see Fig. 
2a). 

Definition 1. A contour C is the ordered set of n edges {Eo, Ej, ... , 
E1}, n23 ,so that Vi € {0...n-1} E_, =E(v;4, vi) and E=E(y;, vi)". 

Point v; for which the edge £;_, (resp. E;) is previous (resp. next), is 
called ith vertex of the contour C, in this case the edges £i- and Æ; are called 
incident. The order of tracing the edges of the contour C = {Ep, E1, ... , En-1} 
by increasing the edge index (i > i + 1) is called straight, the opposite order is 
called reverse. 

Note that the same geometric point in the plane can correspond to several 
vertices of a contour. Such vertices are called the high degree vertices. The 
angle of next (resp. previous) edge E(a,v) (resp. E(v,b)) for a vertex v is 
defined as the directed angle between vectors e; (unit vector of the x-axis) and 
E(v, a) (resp. E(v,b)). 


! Here and furthermore all index ariphmetics is made in modulo n, where 
n is anumber of vertices in currently considered contour. 
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Consider a pair of incident edges £; and E; of a contour C: 0 <i<n- 1. 
A set of points x satisfying the conditions 0 < (visi vix) < Z(Vi+ı Vi Vi-1) and 
0 < |x- vi| < €, is called the inner €-neighborhood of vertex v; and denoted as 
gto). A set of points x satisfying the conditions 0 < <(vi-1 vi x) < Z(Wi-1 Vi Vi+1) 
and 0 < |x — vi| < € is called the outer €-neighborhood of vertex v; and denoted 
as € (vj) (see Fig. 2b). 


From now our setting will be Eucledian plane with usual Cartesian 
system. 


Definition 2. Domain is a finite non-empty polygonal open connected 
set on Eucledian plane, defined with a precision of set of 
measure zero. 


left ep-neighborhood of Eva, b) 


right -neighborhood of Ea, b) 


b) 
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Figure 2 


“Polygonal” means that the domain boundary consists of closed polylines. 
The precision of set of measure zero removes from consideration domains with 
“gaps” and “point holes”. Note that domain can be both single-connected and 
multi-connected set. 


Let dA be the boundary of a domain A, so that the orientations of dA and 
A are coordinated. To obtain bijection between dA and A, we add to the 
traditional boundary representation two important restrictions. Consider 
contour C consisting of n edges. 


Definition 3. Contour C is called the bounding contour of a domain A, 
if it satisfies the following conditions: 


1. CcdaA; 
2. when C is traced in straight order, domain 4 is on the left; 


3. Yie{0...n-1}3 €>0:V xes) xed; 
Vie{0...n-1}Y €>0:4 xee v) xg. 
Definition 4. Polygon is a set of all bounding contours of a domain. 


The outer (resp. inner) contours of A are defined as the contours of its 
polygon with counter-clockwise (resp. clockwise) orientation. Obviously there 
is exactly one outer contour in an arbitrary polygon 


Definition 5. Region is a set of polygons describing a set of non- 
intersecting domains. 
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Figure 3 


Thus we have a bijection between a set of non-intersecting domains (i.e. 
arbitrary polygonal area in the plane) and the set of its bounding contours. The 
examples of the bijection are shown in Fig. 3. 


THE ALGORITHM 


Input: two regions A and B, operation op = {union, intersection, 
difference, symmetrical difference}. 


Output: region R = A op B. 

The algorithm consists of the following steps. 
Step 1. Processing of the edge intersections. 
Step 2. Edge and contour labeling. 
Step 3. Collecting the resulting contours. 
Step 4. Placing the resulting contours in R. 


Fig. 4 shows example regions A and B. According to our definitions, 
region A is a polygon consisting of one outer and one inner contours, and 
region B is a polygon consisting of one outer self-touching contour. 


Figure 4 


Processing of the edge intersections 


Assume we have an algorithm, which reports all intersections between 
non-incident edges of regions A and B. If the edges share common line 
segment, we will treat its endpoints as the intersection points. All new 
intersection points are added as vertices to the input regions. Hence insertion of 
the intersection points does not change domains whose interiors are described 
by A and B, furthermore 4 and B will mean input regions A and B with all their 
intersection points added. Vertices corresponding to intersection points are 
called cross-vertices (Fig. 5). Note that one geometric intersection point 
corresponds to at least two cross-vertices. 

Let v be a cross-vertex of some region. Let E‘(v) (resp. E (v)) denote its 
previous (resp. next) edge. The cross-vertex is processed as follows. 

Two cross-vertex descriptors D*(v) and D™(v) are created, which 


correspond to E*(v) and E (v). Into these descriptors are placed the values of 
their edge angles and flags about if the edges are previous or next to v. 


1. The pointers between v and its descriptors D*(v) and D™(v) are established, 
allowing determining the descriptors for a vertex and vice versa. 


2. D*(v) and D‘(v) are placed into connectivity list L(x) of intersection point 
x. This list contains descriptors of all cross-vertices corresponding to x and is 
sorted by edge angles. Furthermore, L(x) will be used for searching for counter- 
clockwise and clockwise neighbors of an edge. The order of descriptors with 
equal edge angles is not specified. 
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Figure 5 


Fig. 6 shows the connectivity list of point x (marked in Fig. 5). The 
corresponding cross-vertices are as, B, and Bo. 
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Figure 6 


A cross-vertex is considered to belong to a connectivity list when its 
descriptors belong to the list. Since a cross-vertex may belong only to the 
connectivity list of corresponding geometric point, we introduce the notion of a 
cross-vertex connectivity list. 


Edge and contour labeling 
Let C be a bounding contour of a region A or B. Let M be region, which C 
does not belong to. 


Since in step 1 all intersection points were added to the input regions, 
taking into account the definitions 3 — 5 implies the following predicate. 


Predicate 1. After the first step of the algorithm: 
1) Unequal edges of regions Æ and B intersect at their endpoint or not at all; 


2) Every edge of C either coincides with an edge of M (such edges are called 
shared), or lies inside or outside (possibly except its endpoints) M; 


3) IfC does not have cross-vertices, it lies inside or outside M. 


Let E be an edge of a contour C. Then its /abel has one of the following 
values: 


INSIDE — when E (possibly except its endpoints) lies inside M, 
OUTSIDE — when E (possibly except its endpoints) lies outside M, 


SHARED] — when E coincides with an edge from region M with the same 
direction, 


SHARED2 — when E coincides with an edge from region M with the opposite 
direction. 


The Zabel of contour C has one of the following values: 
ISECTED — when C contains at least one cross-vertex, 
INSIDE — when C lies inside M, 

OUTSIDE — when C lies outside M. 


The correctness of edge and contour label definitions immediately follows 
from Predicate 1. 


Here is an algorithm for labeling C and its edges. 
1. C does not contain any cross-vertices. 


If C lies inside M, it is labeled as INSIDE, otherwise — as OUTSIDE. The 
edges of C are not labeled at all. 


2. C contains at least one cross-vertex. 


C is labeled as ISECTED and all its edges are sequentially labeled. Let E;(a, b) 
(i € {0...2-1}, where n is the number of edges in C) be the edge to label. 


2.1 E; does not contain cross-vertices as its endpoints. 


1) i#0. The edge label value is copied from the edge £;_;. 


2) i=0.If a lies inside M, then Æ; is labeled as INSIDE, otherwise Æ; is 

labeled as OUTSIDE. 
2.2 E; contains one or two cross-vertices as its endpoints. 

1) dedge F(c, dd: Fe MAa=can b= d. Eis labeled as SHAREDI. 

2) Aedge F(c, d: Fe Mn a=da b=c. Eis labeled as SHARED2. 

3) Cross-vertices connectivity list(s) does not contain any vertices from 
M. 
Such situation is possible if C is not intersected by M and touches 
itself. The edge is labeled similarly to step 2.1. 

4) Cross-vertices connectivity list(s) contain vertices from M. 
For each such vertex wx the check is performed if £; lies inside 
Z (Wk-1WkķWķ+1). If Æ; is inside one of such “labeling” angles, then it lies 
inside M and therefore is labeled as JVS/DE. Otherwise, Æ; is labeled 
as OUTSIDE (see Fig. 7). 
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Figure 7 
Collecting the resulting contours 


Preliminaries 


One of the essential ideas of the presented algorithm is that bounding 
contours of R (furthermore they are called resulting contours) are collected 
using the values of edge and contour labels instead of their coordinates. Since 
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an edge cannot appear in a region more than once, we are to find the rules for 
edge and contour inclusion into R. 


Fig. 8 shows our example regions A and B after performing labeling step. 


A 
OUTSIDE 
—_- 


INSIDE 
> 


SHARED2 


e 


SHARED]I 


—> 


Figure 8 


The rules for edge inclusion into a resulting contour 


Consider some edge E belonging to either A or B. Let Ọ > 0 be the 


maximum real number such that neither left, nor right @-neighbourhood of 
edge E do not contain any other vertices or edges from A or B. These 
neighbourhoods are called the minimum left and right neighbourhoods of E and 
are denoted respectively as @ (E) and @ (E). The choice of Ọ implies that all 
points from a minimum neighbourhood of an edge are inside or outside A or B 
at the same time. 


Now we show how the rules for edge inclusion into a resulting contour 
follow from those of minimum neighbourhoods. 


1. If @ (E) belongs to R, but @ (E) does not, E will be included into a 
resulting contour with its original direction. 
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2. If @ (E) belongs to R, but @ (E) does not, E will be included into a 
resulting contour with inverse direction. 

3. If both @ (E) and @ (E) do not belong to R, E will not be included into a 
resulting contour. 

4. If both @ (E) and @ (E) belong to R, E will not be included into a 
resulting contour (this would contradict definition 3, clause 4). 


Note that only one edge from a pair of shared edges can be included in 
result. Table 1 shows the inclusion rules for all combinations of edge labels for 
both regions. 


Table 1 


Case 1. 

E belongs to A and its label is 
OUTSIDE. Q` (E) belongs to A, Ọ (E) does 
not belong to any of the input regions. £ will 
be added with the original direction into the 
union, the difference and the symmetrical 
difference of A and B. 

Case 2. 

E belongs to B and its label is 
OUTSIDE. ()(E) belongs to B, 0 (E) does 
not belong to any of the input regions. £ will 
be added with the original direction into the 
union and the symmetrical difference of A 


| and B. 
»- — Case 3. 
- AO . . 
p” da E belongs to A and its label is INSIDE. 
J _ ’‘’ : 
~ F p ( (E) belongs to both A and B, @ (E) 
D i belongs to B. E will be added with the 


original direction into the intersection of A 
and B, and with the inverse direction into the 
symmetrical difference of A and B. 
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Case 4. 

E belongs to B and its label is JNS/DE. 
@ (E) belongs to both A and B, Ọ (E) 
belongs to 4A. E will be added with the 
original direction into the intersection of A 
and B, and with the inverse direction into the 


difference and the symmetrical difference of 
A and B. 


Case 5. 

A pair of edges E € A and F e B with 
labels SHARED2. (p'(E) belongs to A, a 
@ (F) belongs to B. One of the edges will be 
added with the original direction of Æ into the 
difference of A and B. 

Case 6. 

A pair of edges A and B with labels 
SHARED]. The minimum left neighbourhood 
of the edges belongs to both A and B, the 
minimum right neighbourhood of the edges 
does not belong to any of the input regions. 
One of the edges will be added with the 
original direction into the union and the 


| intersection of A and B. 


The jump rules at the intersection points 


Suppose while collecting the resulting edges we just came to a cross- 
vertex v. Where will we go next? Let x be geometric point corresponding to v. 
Consider all cross-vertices corresponding to x. Let m be the number of such 
vertices, m = 2. Then the number of their descriptors is 2m (see Step 1). With a 
help of L(x) we can obtain the set of edges F}, je {0...2m—1} so that Fy = E*(v) 
and the descriptor of Fj,; is the nearest in the clockwise direction to the 


descriptor of F}. 


Let € > 0 be the maximum real number such that no inner or outer €- 
neighbourhood of considered vertices contain vertices or edges different from 


F;. Consider the sectors into which geometric €-neighbourhood of point x is 
divided by edges F;. Such sectors are called the minimum sectors. Apparently, 
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all points of a minimum sector (possibly except its boundaries) are inside or 
outside A or B at the same time. If the points of two neighbour minimum 
sectors belong to R, by definition 3 they should be bounded by the same 
resulting contour. 


We find the desired edge by sequentially testing edges Fj, beginning from 
Fo and increasing j, to satisfy the inclusion rule for desired operation. Then we 
can go along this edge. Note that for the search the geometric information is 
not required at all, because descriptors are sorted in connectivity list by 
corresponding edge angles. 


The collecting algorithm 

Since the previous discussion applies to all operations, the algorithm 
collects resulting contour for all operation in the similar way. In this step we 
sequentially consider each bounding contour of A and B. Let C be the contour 
currently being considered. 

The label of contour or edge X is furthermore denoted as X. Flags. Also 
FORWARD and BACKWARD will denote the original and the inverse 
directions of edges. Here is the algorithm for collecting the resulting contours. 
1. If C.Flags + ISECTED, C is added to R depending on its label and 
operation op (see Table 2). 


Table 2 

op | The rule for inclusion of C into R | Direction 

5 (C.Flags = INSIDE) | FORWARD 

U (C.Flags = OUTSIDE) FORWARD 

= (C. Flags = OUTSIDE) ^ (C € A) FORWARD 
(C.Flags = INSIDE) ^ (C € B) BACKWARD 

x (C.Flags = OUTSIDE) FORWARD 
(C.Flags = INSIDE) BACKWARD 


2. If C.Flags = ISECTED, we need to find in C all edges which resulting 
contours can start from. Here is the algorithm for processing C (n is the number 
of vertices in C). Each edge E has a bit flag E. Mark indicating if the edge has 
already included into one of the resulting contours, initially the flag is set to 
false for all edges. 

for i := 0 to n-1 do 

begin 
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if (EdgeRule(£,, dir) and not £; .Mark) 
begin 
contour r; 
if (dir = FORWARD) 
r := Collect(v;, dir); 
else 
r := Collect(v;,,, dir); 
Include r into a set of resulting contours; 
end; 
end; 
where function EdgeRule(E:edge; var dir:;FORWARD, BACKWARD)): 
boolean is the edge inclusion rule for currently performed operation (see Table 


Table 3 
o | The rule for inclusion of E into R | Direction 
p 
M | (E.Flags = INSIDE) v (E.Flags = SHARED) | FORWARD 
U | (Flags = OUTSIDE) v (E.Flags = SHAREDI) | FORWARD 
— | (EFlags = OUTSIDE) FORWARD 
v (E. Flags = SHARED2)) ^ (E € A) 
((E.Flags = INSIDE) BACKWARD 
v (E. Flags = SHARED2)) ^ (E € B) 
@ | (EFlags = OUTSIDE) FORWARD 
(E.Flags = INSIDE) BACKWARD 


function Collect(v:vertex; dir: FORWARD, BACKWARD)):contour; 
begin 
Create an empty contour r; 
repeat 
Add v tor; 
if (dir = FORWARD) 
E := edge next to v; 
else 
E := edge before v; 
E.Mark := true; 
if (E. Flags = SHARED!) or (E.Flags = SHARED2)) 
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for edge, shared with E, set its Mark to true; 
v := vertex next to v in direction dir; 
if (v is a cross-vertex) then 
Jump(v, dir); 
until (current edge is marked); 
return r; 
end; 


procedure Jump(var v:vertex; var dir: FORWARD, BACKWARD)); 
begin 
if (dir = FORWARD) then 
d := prev(D*(v)); 
else 
d := prev(D (v)); 
{ prev(D) denotes descriptor nearest to D in clockwise direction } 
found := FALSE; 
repeat 
e := edge corresponding to d; 
if (not (e.Mark) and EdgeRule(e, newdir)) then 
begin 
v := vertex corresponding to d; 
if ((e is next to v) and (newdir = FORWARD) or 
(e is previous to v) and (newdir = BACKWARD)) then 
begin 
dir := newdir; 
found := TRUE; 
end; 
end; 
d := prev(d); 
until (found); 
end; 


Fig. 9 shows the resulting contours for different operations on example 
regions A and B. 


The union is composed of one outer and two inner contours, both the 
difference and the intersection are composed of two outer contours, the 
symmetrical difference is composed of three outer contours. 
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Placing the resulting contours in R 


After previous steps, we have obtained a set of bounding contours of R. 
Arranging them into R is very simple because they already have proper 
orientation. 


For each outer contour, a new polygon is created. Inner contours were 
stored in temporary list in previous stage, now they should be placed in proper 
polygons. For finding a proper polygon for an inner contour, it suffices to 
determine the minimum outer contour containing the inner contour. Region R is 
composed of polygons formed in this step. 


AUB AAB 
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THE ANALYSIS OF THE EFFICIENCY 


Now we analyze the efficiency of the presented algorithm. Let A and B be 
input regions with n vertices and z contours in total. 


Let k be the number of new vertices added in step 1. In the worst case k 
= O(n’). Time complexity of step 1 depends on the algorithm for reporting 
edge intersections. The worst-case time complexity of testing all pairs of edges 
for intersection O(n?). There are more efficient algorithms [1. — 3. ] with the 
time complexity O(n log n+ k). Finke and Hinrichs [4. ] presented an algorithm 
which can be extended to compute overlay of two path-connected planar 
subdivisions in O(n log*n+ k + z log n) expected time. 


The time complexity of step 2 is O(n + z p(n)), where p(n) is the time for 
determining if a given point is inside a region consisting of n vertices. Without 
using additional data structures p(n) = O(n). To speed up point location queries 
one can use a search structure created by Seidel’s trapezoidation algorithm, 
which can be constructed in O(n log*n + z log n) time [11. ]. This structure 
allows point location queries in O(/og n) expected time. 


Since in step 3 each edge is used not more than two times, its time 
complexity is apparently O(n + k). Step 4 takes O(z p(n + &)) time. 


Computational expenses of the algorithm can also be decreased by using 
minimum bounding boxes which can be computed in O(n) time. 
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Thus, the presented algorithm works in O(n log*n+ k + z log n) expected 
time and O(n + k) space. 


CONCLUSION 
Let us summarize the main results of this paper. 


e Mathematically non-ambiguous definitions are given for description of 
polygonal regions in the plane. 


e The technique for explicit handling vertices of high degree is developed, 
which avoids degeneracy problems of traditional algorithms. 


e A set of intuitive clear rules for collecting the resulting contours is 
introduced, which allow compact and efficient implementation of the presented 
algorithm. 


The main contribution of the presented paper is a closed set of operations. 


The implementation of the algorithm consists of about 850 lines of C- 
code. In practice, the performance of the algorithm meets theoretical 
requirements very well. 


The first author would like to thank Klamer Schutte for helpful 
discussions on polygon clipping and Ralf Gueting, Raimund Seidel and Jack 
Snoeyink for assistance with materials for the paper. 
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